Study species

Lodgepole pine is a wind-pollinated monecious conifer with an extensive range and well documented local adaptation. While lodgepole pine has four subspecies, this work concerns only Pinus contorta subsp. latifolia. Pollen and seed cone buds differentiate in late summer and early fall, then go dormant. As temperatures warm the following year, buds resume development and strobili “flower”; receptive female strobili exude pollination drops between bracts that capture pollen shed from mature male strobili.

Rangemap of lodgepole pine. After Little E.L. (1971) .

Data

Phenology Data

This project takes advantage of an existing lodgepole pine pollination phenology dataset collected over a decade and a half by government and industry workers in seed orchards. Seed orchards in British Columbia produce large amounts of tree seed for reforestation from parent trees sourced from provenances around the province. To plan for future seed production and orchard establishment and management, seed orchard managers monitored pollination phenology and seed production. Pollination phenology data was collected at the Prince George Tree Improvement Station in British Columbia beginning in 1997 and collection at many other BC tree orchard sites began in 2006 under the Forest Genetics Council’s Operational Tree Improvement Program 0722 (???). Collection continued intermittently through 2012.

Sites

Trees selected from across the British Columbia portion of the lodgepole pine range are grown in seed orchards as part of tree breeding and seed production programs. Between 1997 and 2011, flowering phenology of lodgepole pine was recorded at 7 seed orchard sites in the interior of British Columbia. I contacted Seed Orchard Managers and other forestry professionals across British Columbia in 2012 and received pollination phenology data from C. Walsh, previously at Kalamalka Seed Orchards (now retired), R. Wagner at the Prince George Tree Improvement Station, and J.E. Webber previously at the Glyn Road Research Station (now retired). 4 of the sites are clustered near to one another around Vernon, BC. Sites span about 5 \(^\circ\) of latitude and are at elevations from 466 to 638 m.

Map of seed orchard locations. Seed orchard locations are in red. Weather data gridpoints are in blue. Boxed area in map on left is shown in greater detail on the right.

Provenances

Trees grown at the seed orchard sites were sourced from 7 Seed Planning Units (SPUs). SPUs are organizational units in BC forestry that are based on biogeoclimate regions. Trees with the same SPU provenance are grown together in an orchard at a given site. Genotypes (labelled with a Clone number in the data) in the orchards are represented by multiple ramets.

I obtained pollination phenology records from 17 of the 26 lodgepole pine seed orchards in British Columbia. Orchards include the offspring of trees from 5 of the 13 BC seed planning zones (SPZs) [SPZ REFERENCE OR FIGURE?] grown at at 7 sites across BC. SPZs are biogeoclimatic and political units used for seed planting purposes by British Columbia. SPZs are divided into elevation bands called Seed Planning Units (SPUs), which form this project’s 7 provenances .

Thirteen out of 17 orchards in my data set are first generation orchards and should faithfully represent their provenances. These first generation orchards represent 6 provenances at 5 sites. Second generation orchards have been selectively crossed and this may skew the mean or variance of phenology for a provenance if pollination phenology varies by provenance.

Most provenances are represented at 2 to 3 sites and have at least three years of data at a given site spanning 1997-2012 (Figure ). The Prince George Tree Improvement Station (PGTIS) provides a continuous 15-year record of its three orchards’ phenology. Most sites only observe one ramet per clone, PGTIS typically observed 2-4 ramets per clone.

Contingency table of years of data for Seed Planning Units (rows) and Seed Orchard Sites (columns). Seed Planning Zones, used as provenances in this project, are usually represented at multiple years and multiple sites. There is particularly good representation at PGTIS.

Map of Seed Planning Units (SPUs). Seed planning units are biogeoclimatic and political units used for seed planting purposes by British Columbia. Seed planning units form this project’s provenances. High, Low, and Mid refer to elevational bands. Data is also available for East Kootenay Low, but will likely not be included in any analysis as it includes only one year at one site.

[table with Site Columns and SPU rows with years of data as values] [table of number of trees/clones in a given year for a given site/spu combo]

Phenology scoring protocol

Protocol C in (Woods, Stoehr, and Webber 1996) was used as the basis for collecting pollen shed and cone receptivity data, though operational constraints led to some modifications. Workers monitored seed orchards for the beginning of pollen shed and cone receptivity. ~15 clones, usually represented by 2 trees each, were selected for specific observations. When the active period seemed to be starting, workers went into the orchard every few days to make observations on the selected trees.

Stages of pollen and seed cone development are described by [Owens & Molder 1984 and updated in Owens2006] and were used as a general guide for determining the phenological state of pollen and seed cones. Pollen cones are flowering when tapping causes pollen to be released and seed cones are flowering when there are gaps between the bract-scale complexes. Pollination drops are also produced during flowering, though they recede midday if pollen is not present (???). q8

“Flowering” states were recorded for both pollen shed and cone receptivity at the level of each tree. Protocol C recommends marking the dates when 20% of the cones on a tree have begun flowering and when 80% of the cones on a tree have finished flowering. Operationally, there was some subjectivity and tree-level states for each cone type should be interpreted as “starting flowering” and “finished flowering.” There is some subjectivity here.

[DESCRIBE SURVEY PERIODS]

# Index by Site, Provenance, and Year
dat$spyind <- group_indices(dat, Site, Year, Provenance)

surveydf <- select(dat, spyind, Year, Site, Provenance, DoY) %>%
  distinct()

ggplot(surveydf, aes(x=DoY, y=as.factor(spyind), color=Provenance, shape=Site)) +
  geom_point() +
  geom_line() +
 facet_grid(rows=vars(Site,Year), scales="free_y") +
  scale_color_viridis_d(option="C") +
  scale_shape_manual(values=c(1:7)) +
  theme_bw(base_size = 15) +
  theme(strip.text.y = element_text(angle = 0), axis.text.y = element_blank(), axis.ticks = element_blank(), axis.title.y = element_blank(), panel.border=element_blank(), panel.grid.minor.y=element_blank(), panel.grid.major.y = element_blank()) +
  ggtitle("Observation Dates")

Observations were not made every day and survey periods varied in length. At Prince George, not all trees have complete phenological records, e.g if a tree is not flowering on the first day of observation, the start date is unknown.

Harmonization and cleaning

[DESCRIBE DATA TRANSCRIPTION AND CLEANING].

There were some differences in how data was recorded at the Prince George Tree Improvement Station versus the other sites. At Prince George, trees were marked as flowering or not flowering on each day of observation. At other sites, only the first day observed flowering and the first day observed finished flowering were recorded. I inferred that trees were not yet flowering on observation days prior to their first recorded flowering date. I cleaned and harmonized the data for analysis in a single model using R scripts provided so that, for pollen cones and seed cones, stage 1 = not yet flowering, stage 2 = flowering, and stage 3 finished flowering.

There are three phenological stages of interest each for a tree’s male and female cones

• 1. The cones have not yet flowered
• 2. The cones are flowering
• 3. The cones have finished flowering

Phenophases in the field were recorded using different symbol sets and resolutions. I assigned each symbol to one of the phenophases above. Trees that did not produce cones are assigned phenological stage 0.

In total, there were 21909 usable observations of phenological state.

Some observations at the Prince George Tree Improvement Station Site were redundant and were dropped for modeling. For example, if a tree’s pollen cones were recorded as “finished flowering” for 3 days, only the first observation was kept - phenological states are not reversible. This improved the speed of model fitting.

Phenophase	Symbols	Male.Cones	Female.Cones
0	0	none produced	none produced
1	1, 2.5, -	not yet shedding	not yet receptive
2	3, 3.5, 4, 4.5, 5, pollenshed20, receptive20	shedding	receptive
3	-, receptive80, pollenshed80	finished shedding	no longer receptive

[Figure showing phenology data somehow]

Weather data

Daily weather data at seed orchard sites was extracted from PNWNAmet, a daily gridded meteorological dataset at 1/16 \(^\circ\) over northwest North America (Consortium 2014). The closest point in the PNWNAmet grid was used for each station . Mean daily temperature was calculated as the average of the minimum and maximum daily temperatures provided by PNWNAmet.

As the PNWNAmet gridpoints do not exactly align with the sites, a correction for elevation was applied using the dry adiabatic lapse rate.

\[t_{site} = t_{gridpoint} + \left(\frac{10^{\circ}\mathrm{C}}{1000 \mathrm{m}}\right)\left(elevation_{site} - elevation_{gridpoint}\right)\]

FIGURE SUMMARIZING RAW MEAN TEMPERATURES?

Figure X shows the corrected and uncorrected temperatures with temperatures collected by seed orchard workers at the sites using automated data loggers. In the figure, PNWNAmet temperatures have been translated to growing degree days with a 5\(^{\circ}\)C threshold because only GDD calculated from the logged data was available from the sites.

Forcing units

Flowering is a developmental process that speeds up and slows down according to the current temperature. Forcing units describe the relative effect of temperature on development. Observable phenological events occur only after a certain amount of forcing has accumulated. The mean daily temperature \(t\) was mapped to forcing units with a function that describes the relationship between temperature and development rate.

\[R(t) = \frac{1}{1+ e^{-0.185*t-18.4}}\] This equation for forcing units is based on experimental work in temperate forest tree species (Sarvas 1972, @hanninenModellingBudDormancy1990) and is verified to perform better than growing degree day models by Chuine (Chuine, Cour, and Rousseau 1999).

Accumulated forcing units on day \(d\) (\(f(d)\)) are the sum of the relative temperature effect from January 1 (ordinal date \(1\)) to day \(d\).

\[f(d) = \sum_{i=1}^d R_d(x)\]

Cooler sites, like Prince George, accumulate forcing slower than warmer sites like Kalamalka.

[ACCUMULATED FORCING UNITS FIGURE]

Phenology Modeling

I modeled phenological states as a function of accumulated forcing units with a Bayesian multilevel ordinal logistic model fit in Stan with the No-U-Turn Sampler (Team 2018). The model was fit separately to data for male strobili and female strobili. A logistic cumulative link function accounts for the ordering of phenological states and relates phenological states to a linear model. This type of model makes no assumption about the distance between phenological states. The model parameters describe a probability function for each state. Cutpoints \(\kappa\) separate each state while the slope parameter of the linear model influences the steepness of the curves. Transitions between phenological states happen rapidly when the slope is large and slowly when it is small.

The likelihood of being in state \(s \in \{1,2,3\}\) is

\[\mathrm{Pr}(S=s)=\text{OrderedLogistic}(s~|~\eta,\kappa) = \left\{ \begin{array}{ll} 1 - \text{logistic}(\eta - \kappa_1) & \text{if } s = 1 \\[4pt] \text{logistic}(\eta - \kappa_1) - \text{logistic}(\eta - \kappa_2) & \text{if } s=2 \\[4pt] \text{logistic}(\eta - \kappa_2) - 0 & \text{if } s = 3 \end{array} \right.\]

where the cutpoints \(\kappa_s < \kappa_{s+1}\) and \(\eta\) is a linear model.

[CORRECT THIS TO MATCH STAN CODE] \[\begin{array}{rlr} S_i & \sim \mathrm{OrderedLogistic}(\eta,\kappa) & \text{probability of data}\\ \eta_i & = (\beta + \beta_{site} +\beta_{prov} + \beta_{clone} + \beta_{year}) f_i & \text{linear model}\\ \kappa_k & \sim \mathrm{gamma}(7.5,1) & \text{fixed priors}\\ \beta & \sim \mathrm{exponential}(2) & \text{}\\ \beta_{site} & \sim \mathrm{normal}(0,\sigma_{site}) & \text{adaptive priors}\\ \beta_{prov} & \sim \mathrm{normal}(0,\sigma_{prov}) & \text{}\\ \beta_{year} & \sim \mathrm{normal}(0, \sigma_{year}) \\ \beta_{clone} & \sim \mathrm{normal}(0, \sigma_{clone}) \\ \sigma_{site} & \sim \mathrm{exponential}(3) & \text{hyperpriors}\\ \sigma_{prov} & \sim \mathrm{exponential}(3) \\ \sigma_{year} & \sim \mathrm{exponential}(3) \\ \sigma_{clone} & \sim \mathrm{exponential}(3) \\ \end{array}\]

The linear model \(\eta\) includes a population mean slope \(\beta\) with deviations from that for site, provenance, clone and year. Priors on \(\kappa\)s and \(\beta\) force the correct sign (positive) and order of magnitude. Priors on the site, provenance, clone, and year effects are adaptive and effects are estimated from the data with partial pooling.

Low variation for provenances and years made sampling inefficient, so I used a non-centered parameterization in the Stan model code for these clusters. The model was run for [X] iterations of which [Y] were warmup [See code].

Model parameters

fstart and fend

Posterior predictive checks

[See p. 81 in lab notebook for figure idea]

Average Predictive Comparison

Provenance differences

[FIGURE SHOWING THE PROPORTION OF TREES THAT HAVE REACHED STAGE 2 AND STAGE 3 AT EACH SITE/PROV COMBINATION over forcing units]

1. How strong are the genetic vs. environmental effects on pollination phenology?

To answer this question I will compare the parameters of models fit by provenance. If different provenances have significantly different model parameters - that is, have different heat sum requirements for pollen shed and cone receptivity by provenance, then the rest of my project will be limited to the provenances I have data from. If not, I can use the same models across the entire range when predicting phenological events.

The budburst and spring pollination phenology of temperate tree species is determined by temperature, and there may be some regional variation in the temperatures required for pollen shed and cone receptivity; in lodgepole pine, the threshold for shoot elongation was 5.1 \(^\circ\mathrm{C}\) for coastal and interior provenances, but 6.5 \(^\circ\mathrm{C}\) for northern provenances (Chuine, Aitken, and Ying 2001).

However, seedling spring vegetative phenology data from four growth chamber temperature and moisture treatments in the AdapTree project suggests that budburst shows substantial clinal variation (> 10% of phenotypic variance explained) only in the coldest treatment tested (MAT of 1 \(^\circ\mathrm{C}\)), and much colder than the seed orchard locations (Figure , Figure , Liepe (2014) (submitted)). Analysis of phenotypic data from an outdoor common garden in Vancouver (personal commmunication, Ian Maclachlan, AdapTree Graduate Student, January, 1 2015) also shows that very little variation exists among provenances for spring phenology, regardless of whether seedlings originate from wild-stand or selectively bred seed orchard seedlots.

Variance in budbreak and budset for lodgepole pine explained by provenance climate. Seedlings from 281 locations across lodgepole’s BC range were grown in growth chambers under four treatments: approximated seasonal and daily variation in temperature for geographic locations with a MAT of 1, 6, and 11 \(^\circ\mathrm{C}\), with the warmest treatment having a dry and wet version. Pearson’s correlation coefficient between phenotypic traits in lodgepole pine and provenance climates are presented. The vertical line represents the critical \(r^2\) value after an adjustment for multiple inferences. From Liepe (2014) (submitted).

2. How does faster spring warm-up like we expect under climate change affect within population mating success and frequency of outcrossing?

I expect that faster spring warm-up will condense as well as advance growth initiation dates and also shorten pollination phenological periods by speeding development prior to and during the phenological period. This should further synchronize populations in similar climates and act as a barrier to breeding with populations in more different climates. On the other hand, it could favor outcrossing over selfing by causing protandry - pollen shed in advance of cone receptivity. Based on anecdotal reports from seed orchards managers, I expect protandry (and thus outcrossing) to increase under climate change.

I will examine the variability of the start date across the range (\(t_1\)) to see if it shrinks in warmer years and as the climate changes. I will also calculate overlap between within population cone receptivity and pollen shed under current conditions and with climate warming and see if the heat sum requirements for cone receptivity and pollen shed are different.

Model summary

# ffit.stan <- readRDS("../FEMALE_slopes_scaled.rds")
# mfit.stan <- readRDS("../MALE_slopes_scaled.rds")
# 
# fsum <- rstan::summary(ffit.stan)$summary
# fsum <- as.data.frame(fsum)
# fsum$par <- rownames(fsum)
# fpost <- as.data.frame(ffit.stan)
# fpost_re <- extract.samples(ffit.stan)
# fparam_names <- colnames(fpost)
# 
# flp <- log_posterior(ffit.stan)
# fnp <- nuts_params(ffit.stan) #nuts params
# 
# msum <- rstan::summary(mfit.stan)$summary
# msum <- as.data.frame(msum)
# msum$par <- rownames(msum)
# mpost <- as.data.frame(mfit.stan)
# mpost_re <- extract.samples(mfit.stan)
# mparam_names <- colnames(mpost)
# 
# mlp <- log_posterior(mfit.stan)
# mnp <- nuts_params(mfit.stan) #nuts params
# 
# #mf
# fpost$sex <- 0
# mpost$sex <- 1
# 
# 
# post <- full_join(mpost, fpost) %>%
#     select(-`b_clone[260]`)

# kable(
# fsum[str_detect(fsum$par, "b_clone", negate=TRUE) , ], digits=3)

# kable(
# msum[str_detect(msum$par, "b_clone", negate=TRUE) , ], digits=3)

# compare_fm <- function(femplot, mplot, nrow = 2, ...) {
#     bayesplot_grid(
#         femplot, mplot,
#         grid_args = list(nrow = nrow),
#         subtitles = c("Female",
#                       "Male"),
#         ...
#     )
# }
# 
# compare_fm_wide <- function(femplot, mplot, ncol = 2, ...) {
#     bayesplot_grid(
#         femplot, mplot,
#         grid_args = list(ncol = ncol),
#         subtitles = c("Female",
#                       "Male"),
#         ...
#     )
# }

# ggplot(post, aes(x=beta, fill=as.factor(sex))) +
#     geom_density(alpha=0.5) + 
#     ggtitle("male and female speed of transition estimates")
# 
# fint_beta <- mcmc_areas(fpost, regex_pars = "beta") + ggtitle("beta")
# mint_beta <- mcmc_areas(mpost, regex_pars = "beta")
# 
# compare_fm(fint_beta, mint_beta, xlim=c(0.05, 0.10))
# 
# fkap <- mcmc_areas(fpost, regex_pars = "kappa")
# mkap <- mcmc_areas(mpost, regex_pars = "kappa")
# 
# compare_fm_wide(fkap, mkap, xlim=c(18, 29))
# 
# color_scheme_set("blue")
# fint_site <- mcmc_intervals(fpost, regex_pars = c("site", "prov"))
# mint_site <- mcmc_intervals(mpost, regex_pars = c("site", "prov"))
# 
# compare_fm_wide(fint_site, mint_site) 
# 
# fa_sigma <- mcmc_areas(fpost, regex_pars= "sigma") + ggtitle("variance")
# ma_sigma <- mcmc_areas(mpost, regex_pars= "sigma")
# 
# compare_fm_wide(fa_sigma, ma_sigma)
# 
# fa_y <- mcmc_areas(fpost, regex_pars = "year") + ggtitle("year")
# ma_y <- mcmc_areas(mpost, regex_pars = "year")
# 
# compare_fm_wide(fa_y, ma_y)

# fint_c <- mcmc_intervals(fpost, regex_pars = c("clone"), point_est="none") + ggtitle("Clone effects")
# mint_c <- mcmc_intervals(mpost, regex_pars = c("clone"), point_est="none")
# 
# compare_fm_wide(fint_c, mint_c)
# ```
# 
# ### Sex diff
# 
# ```{r half transitions - no effects}
# fpost$f1 <- fpost$`kappa[1]`/fpost$beta
# fpost$f2 <- fpost$`kappa[2]`/fpost$beta
# mpost$f1 <- mpost$`kappa[1]`/mpost$beta
# mpost$f2 <- mpost$`kappa[2]`/mpost$beta
# 
# ff <- mcmc_areas(fpost, pars = c("f1", "f2")) + ggtitle("half transitions")
# mf <- mcmc_areas(mpost, pars = c("f1", "f2")) + ggtitle("half transitions")
# 
# compare_fm(ff, mf)

Limitations of the dataset

Benefits of this dataset include the large number of trees, long time series, multiple provenances grown at multiple sites giving a semi-common garden design, and the inclusion of clones. Limitations include interval and end censoring, especially at Prince George, irregular scoring systems, subjective scoring, non-random clone selection, and selective breeding.

Grafting bias?

Problems from breeding

There are two types of orchards: first generation and advanced. Each orchard contains trees that are the descendants of seeds and scions collected from mother trees in one particular provenance. First generation orchards have not undergone selective breeding but may be subject to selection through the process of choosing mother trees (called “plus” trees by tree breeders), testing and selection of mother tree offspring, or culling of unfavorable trees in the orchard. First generation orchard trees are clones of mother trees and their offspring grafted onto rootstock. There are three types of first generation orchards subject to different levels of selection: 1.0, 1.5, and 1.75 generation orchards. All orchards are subject to selection based on mother tree selection. In 1.5 generation orchards, superior offspring of mother trees are selected for the orchard or poor trees are culled from the orchard. This selection may be strong; more than 50% of offspring can be rejected in this process due to poor growth form (Ukrainetz 2011). In 1.75 generation orchards, data from 10 year tests of mother tree offspring are used to select 40 of the best genotypes for a given provenance. Offspring from controlled crosses between trees in first generation orchards are tested in field trials and the best trees are then cloned and grafted onto root stock in advanced generation orchards.

orchards_in_data <- read.csv("~/Documents/research_phd/data/PhenologyAndPollenCounts/data_formatted_and_derived/derived_phenophase.csv", header=TRUE, stringsAsFactors = FALSE) %>% 
    select(Orchard, Year, ) %>%
    group_by(Orchard) %>%
    dplyr::summarise(Years = n_distinct(Year))
orchard_gen <- read.csv("~/Documents/research_phd/data/OrchardInfo/OrchardGen.csv", header=TRUE, stringsAsFactors = FALSE)

orchard_meta <- dplyr::left_join(orchards_in_data, orchard_gen) %>%
    unique()

## Joining, by = "Orchard"

knitr::kable(orchard_meta)

selection only from a portion of the range

Provenances and sites included in my data exclude the coldest and wettest parts of the range. Figure shows the 1961-1990 climate normal mean annual temperature (MAT) and mean annual precipitation (MAP) for data from all lodgepole pine provenances in BC, provenances included in my data set, and sites where provenances were grown. Orchard sites are generally much warmer and drier than many locations where lodgepole pine grow. I may have difficulty extrapolating into the coolest and wettest parts of the lodgepole pine range, depending on actual yearly conditions at sites for which I have data. However, this may be advantageous when it comes to predicting phenology under climate change.

Provenance (SPU) climates of data included in this proposal, provenances not included in this proposal, and sites where provenances were grown. MAT is mean annual temperature and MAP is mean annual precipitation.

no chilling

Previous models for lodgepole pine have fitted spring growth phenology models without considering chilling (Chuine, Aitken, and Ying 2001) and, in general, budburst models for boreal tree species have not required chilling to give accurate predictions (Linkosalo 2000). For lodgepole pine growing in natural environments, chilling is always met and should continue to be met under the next century of climate change; in the AdapTree project, lodgepole pine seedlings grown for three to four weeks at 4 \(^\circ\mathrm{C}\) met their chilling requirements. Thus, I expect to be able to disregard the state of chilling and fit only Equations and ). Once the model is parameterized I will be able to calculate the timing of pollen shed and cone receptivity across the lodgepole pine species range using location specific climate data from seed orchards and ClimateWNA (Wang et al. 2012). I will validate the model by comparing predictions to phenological measurements from the National Phenology Network.